{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import graphviz\n",
    "import lingam\n",
    "from lingam.utils import make_dot\n",
    "import logging\n",
    "import dowhy\n",
    "from dowhy import CausalModel\n",
    "from semopy import Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Path analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_excel(r\"C:\\df_analysis_2.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the SEM model\n",
    "model_desc = \"\"\"\n",
    "Distrust ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "Exploitation ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "Indifference ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "Loneliness ~ NDVI + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "\"\"\"\n",
    "\n",
    "\n",
    "# Create and fit the model\n",
    "model = Model(model_desc)\n",
    "results = model.fit(df)\n",
    "\n",
    "# Inspect the model's results\n",
    "params = model.inspect()\n",
    "\n",
    "# Display the results\n",
    "print(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from semopy import Model\n",
    "import graphviz\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "\n",
    "# Define and fit the SEM models\n",
    "model_descriptions = [\n",
    "    \"\"\"\n",
    "    Distrust ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Exploitation ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Indifference ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Loneliness ~ EVI + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    \"\"\",\n",
    "    \"\"\"\n",
    "    Distrust ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Exploitation ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Indifference ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Loneliness ~ NDVI + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    \"\"\",\n",
    "    \"\"\"\n",
    "    Distrust ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Exploitation ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Indifference ~ Loneliness + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    Loneliness ~ LAI + income_2 + education + homeownership + age + sex + marital + religion + ideology\n",
    "    \"\"\"\n",
    "]\n",
    "\n",
    "focus_vars_list = [\n",
    "    [\"EVI\", \"Loneliness\", \"Distrust\", \"Exploitation\", \"Indifference\"],\n",
    "    [\"NDVI\", \"Loneliness\", \"Distrust\", \"Exploitation\", \"Indifference\"],\n",
    "    [\"LAI\", \"Loneliness\", \"Distrust\", \"Exploitation\", \"Indifference\"]\n",
    "]\n",
    "\n",
    "# Function to visualize the SEM model with parameter estimates\n",
    "def visualize_sem_with_focused_params(model, params, focus_vars, path=\"sem_graph\", width=\"80\", height=\"180\", fontsize=\"18\"):\n",
    "    dot = graphviz.Digraph(comment='Focused SEM Model')\n",
    "    dot.attr(size=f\"{width},{height}\", ratio='auto')\n",
    "\n",
    "    # Define the set of variables to focus on\n",
    "    focus_set = set(focus_vars)\n",
    "\n",
    "    # Adding nodes for the focused variables\n",
    "    for var in focus_set:\n",
    "        dot.node(var, var, shape='ellipse')\n",
    "\n",
    "    # Adding edges with parameter estimates for the focused variables\n",
    "    for index, row in params.iterrows():\n",
    "        if row['op'] == '~' and {row['lval'], row['rval']}.issubset(focus_set):\n",
    "            estimate = round(row['Estimate'], 3)\n",
    "            p_value = round(row['p-value'], 3)\n",
    "            if p_value < 0.01:\n",
    "                significance = '***'\n",
    "            elif p_value < 0.05:\n",
    "                significance = '**'\n",
    "            elif p_value < 0.1:\n",
    "                significance = '*'\n",
    "            else:\n",
    "                significance = ''\n",
    "            label = f'est: {estimate}, p: {p_value} {significance}'\n",
    "            dot.edge(row['rval'], row['lval'], label=label, fontsize=fontsize)\n",
    "\n",
    "    # Save and render the graph\n",
    "    dot.render(path, format='png', cleanup=True)\n",
    "\n",
    "# Fit and visualize each model\n",
    "model_paths = []\n",
    "for i, model_desc in enumerate(model_descriptions):\n",
    "    model = Model(model_desc)\n",
    "    results = model.fit(df)\n",
    "    params = model.inspect()\n",
    "    path = f\"focused_sem_example_{i+1}\"\n",
    "    visualize_sem_with_focused_params(model, params, focus_vars_list[i], path, width=\"80\", height=\"180\", fontsize=\"18\")\n",
    "    model_paths.append(f\"{path}.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "# Display the visualizations in a 1x3 subplot layout\n",
    "fig, axes = plt.subplots(1, 3, figsize=(24, 10))\n",
    "\n",
    "# Titles for each subplot\n",
    "titles = [\"(a) EVI and Loneliness\", \"(b) NDVI and Loneliness\"]\n",
    "\n",
    "for i, ax in enumerate(axes):\n",
    "    img = mpimg.imread(model_paths[i])\n",
    "    ax.imshow(img)\n",
    "    ax.set_title(titles[i], fontsize=20)\n",
    "    ax.axis('off')\n",
    "\n",
    "plt.tight_layout(pad=5.0)\n",
    "\n",
    "# Save the figure as a TIFF file\n",
    "plt.savefig(r\"C:\\figure_pathanalysis.tiff\", format='tiff', dpi=300)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
